

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z analyst specif


%% estimated mean utilities
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);


%% counterfactuals
load('counterfactualFR_replication.mat','srng_cfactFR')
rng(srng_cfactUS)
NS = 2e+3;
YY4 = YY(YY(:,6)==0,1:5); % 4-mpow cw-intervention strategy profiles
timeFR=0;

% FR never intervenes
VnoFR_R = UStar(YY(:,4)==0,1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VnoFR_US = UStar(YY(:,4)==0,2,1:N) + randn(size(YY4,1),NS,N);
VnoFR_UK = UStar(YY(:,4)==0,3,1:N) + randn(size(YY4,1),NS,N);
VnoFR_France = UStar(YY(:,4)==0,4,1:N) + randn(size(YY4,1),NS,N);
VnoFR_Russia = UStar(YY(:,4)==0,5,1:N) + randn(size(YY4,1),NS,N);
VnoFR_China = UStar(YY(:,4)==0,6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
noFRintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameFR',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VnoFR_US(2:end,b,n),VnoFR_UK(2:end,b,n),VnoFR_Russia(2:end,b,n),VnoFR_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VnoFR_R(2:end,b,n)'*EQ(:,m);
        end
        noFRintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeFR=timeFR+toc;
end
[noFRDMES,noFRintervEq,noFReqPayoffs] = createDMES(noFRintervEq,VnoFR_R,VnoFR_US,VnoFR_UK,VnoFR_France,VnoFR_Russia,VnoFR_China,[YY4(:,1:3),zeros(size(YY4,1),1),YY4(:,4:5)],specif);

% FR always intervenes
VallFR_R = UStar([1;find(YY(:,4)==1)],1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VallFR_US = UStar([1;find(YY(:,4)==1)],2,1:N) + randn(size(YY4,1),NS,N);
VallFR_UK = UStar([1;find(YY(:,4)==1)],3,1:N) + randn(size(YY4,1),NS,N);
VallFR_France = UStar([1;find(YY(:,4)==1)],4,1:N) + randn(size(YY4,1),NS,N);
VallFR_Russia = UStar([1;find(YY(:,4)==1)],5,1:N) + randn(size(YY4,1),NS,N);
VallFR_China = UStar([1;find(YY(:,4)==1)],6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
allFRintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameFR',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VallFR_US(2:end,b,n),VallFR_UK(2:end,b,n),VallFR_Russia(2:end,b,n),VallFR_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VallFR_R(2:end,b,n)'*EQ(:,m);
        end
        allFRintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeFR=timeFR+toc;
end
[allFRDMES,allFRintervEq,allFReqPayoffs] = createDMES(allFRintervEq,VallFR_R,VallFR_US,VallFR_UK,VallFR_France,VallFR_Russia,VallFR_China,[YY4(:,1:3),ones(size(YY4,1),1),YY4(:,4:5)],specif);


%%
clearvars -except N NS srng_cfactFR YY4 timeFR VnoFR_R VnoFR_US VnoFR_UK VnoFR_France VnoFR_Russia VnoFR_China ...
     VallFR_R VallFR_US VallFR_UK VallFR_France VallFR_Russia VallFR_China noFRintervEq noFRDMES noFReqPayoffs allFRintervEq allFRDMES allFReqPayoffs
save('counterfactualFR.mat')



